In [11]:
import numpy as np
import xarray as xr
import pandas as pd
import scipy
import sys
import warnings
warnings.filterwarnings('ignore')
import geopandas as gpd
from shapely.geometry import mapping
import cartopy.crs as ccrs
import cartopy.feature
import matplotlib.pyplot as plt
import time as timer
start_all = timer.time()
In [12]:
YR = 1994
dataf ="/Volumes/ESA_F4R/era/"
#datao ="/Volumes/ESA_F4R/ed_prepare/"
datao ="/Users/ellendyer/Desktop/"
datap ="/Users/ellendyer/Library/Mobile Documents/com~apple~CloudDocs/1SHARED_WORK/Work/3_ESA_GRANT/MODEL/plots/era/"
datas ="/Users/ellendyer/Library/Mobile Documents/com~apple~CloudDocs/1SHARED_WORK/Work/3_ESA_GRANT/MODEL/Shapefiles/"
shp_cod = gpd.read_file(datas+"geoBoundaries-COD-ADM0.shp")
Read in pre-processed files that have the following conversions:¶
Read in ERA5 data on pressure levels (hourly timesteps in fortnightly files)
- fortnightly files currently run from 1994-2024
- resampled to monthly MS timestep
- shum multiplied by 1000 to convert from kg/kg --> g/kg
- pressure levels are divided by 100 to convert from Pa to hPa (only for fortnightly files)
- sort data by descending pressure levels (only for fortnightly files)
Input file units:
- plev - pa
- q - kg/kg
- u - m/s
- v - m/s
Read in ERA5 land data (hourly in monthly files)
- selecting hour 23 (0-23) of Prec and Evap because of how ERA5 Land variables are accumulated (https://confluence.ecmwf.int/pages/viewpage.action?pageId=197702790 - https://confluence.ecmwf.int/display/CKB/ERA5-Land%3A+data+documentation#ERA5Land:datadocumentation-accumulationsAccumulations)
- prec is multiplied by 1000 to convert from m to mm
- evap is multiplied by -1000 to convert from m to mm and upward fluxes in land model are considered negative
- Prec, Evap, and Psfc are then resampled to MS monthly and also interpolated to coarser pressure level grid
Input file units:
- tp - m
- e - m (-)
- sp - pa
Merging all input datasets into one dataset for recyling code called ds
- close both input datasets
- sort everything so latitude is south to north
- transpose dimensions so they run (lon,lat,level,time) as in recycling code
- save input ds to file
Integrate zonal and meridional moisture flux
- To do: do this with hourly data to test difference in rho
- Must check if there are any nans in input arrays - there can be none because we are using nans as an indicator in the modified definitions
In [13]:
ds = xr.open_dataset(datao+"full_reg_erads_"+str(YR)+".nc")
ds = ds.transpose("time","level","lat","lon",missing_dims='ignore')
#Running the code with an irregular boundary means that there will be nans once the data is clipped
#For this to work properly it is important that there are no nans in the input data which would
#actually be missing data where recycling will be calculated
print(np.isnan(ds['Evap'].values).any())
print(ds['Evap'].isnull().count().values)
#print('Number of nans in Evap: ', ds['Evap'].where(np.isnan(ds['Evap'])==True,drop=True).count().values)
#print('Number of nans in Prec: ', ds['Prec'].where(np.isnan(ds['Prec'])==True,drop=True).count().values)
#print('Number of nans in Psfc: ', ds['Psfc'].where(np.isnan(ds['Psfc'])==True,drop=True).count().values)
#print('Number of nans in Uwnd: ', ds['Uwnd'].where(np.isnan(ds['Uwnd'])==True,drop=True).count().values)
#print('Number of nans in Vwnd: ', ds['Vwnd'].where(np.isnan(ds['Vwnd'])==True,drop=True).count().values)
#print('Number of nans in Shum: ', ds['Shum'].where(np.isnan(ds['Shum'])==True,drop=True).count().values)
#Set evaporation outside local box
mask_lon = (ds['Evap'].lon >= 18) & (ds['Evap'].lon <= 30)
mask_lat = (ds['Evap'].lat >= -9) & (ds['Evap'].lat <= 3)
ds['Evap_local'] = ds['Evap'].where(mask_lon & mask_lat, 0.0)
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cartopy.feature.OCEAN)
#shp_cod.plot(ax=ax, edgecolor='yellow',facecolor='none',lw=2,zorder=2,linestyle='--')
ds['Evap'].mean('time').plot(ax=ax,transform=ccrs.PlateCarree(),
add_colorbar=True,
#vmin=0,vmax=2,
alpha=1,
cmap=plt.cm.viridis,
extend="both")
ds['Evap_local'].mean('time').plot(ax=ax,transform=ccrs.PlateCarree(),
add_colorbar=True,
#vmin=0,vmax=2,
alpha=0.3,
cmap=plt.cm.RdBu,
extend="both")
ax.set_extent([11, 34, -14, 8])
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='black', alpha=0.5, linestyle='dotted')
gl.top_labels = False
gl.right_labels = False
#ax.set_title()
#plt.savefig()
plt.show()
plt.clf()
plt.close()
print(np.isnan(ds['Evap_local'].values).any())
print(ds['Evap_local'].isnull().count().values)
ds = ds.transpose("lon","lat","level","time",missing_dims='ignore')
ds = ds.sel(lon=slice(15,33),lat=slice(-12,6))
True 88740
False 88740
In [14]:
#Prepping datasets near surface for recycling
import bulk_recycling_model.numerical_integration
# Integrate 10^-3 Shum Uwnd dp
# Because the integration limits are from high pressure to low pressure, we need to invert the sign.
integrand = -1 * 1e-3 * ds["Shum"] * ds["Uwnd"]
Fx = bulk_recycling_model.numerical_integration.integrate_no_extrapolation(integrand, ds["Psfc"])
# Units: mb x m/s
# Integrate 10^-3 Shum Vwnd dp
# Because the integration limits are from high pressure to low pressure, we need to invert the sign.
integrand = -1 * 1e-3 * ds["Shum"] * ds["Vwnd"]
Fy = bulk_recycling_model.numerical_integration.integrate_no_extrapolation(integrand, ds["Psfc"])
# Units: mb x m/s
Prepare scaled data for recycling code
- Evaporation and moisture fluxes
In [15]:
# Prepare and scale the data
from bulk_recycling_model import preprocess
from bulk_recycling_model import ED_preprocess
from bulk_recycling_model.axis import Axis
from bulk_recycling_model.scaling import Scaling, UnitSystem
# degrees
L = ds.coords["lon"].max().item() - ds.coords["lon"].min().item()
# convert to meters
L = L * 111e3 * np.cos(np.deg2rad(ds.coords["lat"].mean().item()))
dx = L / ds.sizes["lon"]
# lon axis
lon_axis = Axis(
ds.coords["lon"].min().item(),
ds.coords["lon"].diff("lon").mean().item(),
ds.sizes["lon"],
)
# degrees
H = ds.coords["lat"].values[-1] - ds.coords["lat"].values[0]
# convert to meters
H = H * 111e3
dy = H / ds.sizes["lat"]
# lat axis
lat_axis = Axis(
ds.coords["lat"].min().item(),
ds.coords["lat"].diff("lat").mean().item(),
ds.sizes["lat"],
)
print(f"{L = :.2e} m")
print(f"{dx = :.2e} m")
print(f"{H = :.2e} m")
print(f"{dy = :.2e} m")
# make a scaling object to convert between unit systems
scaling = Scaling(H)
dx = scaling.distance.convert(dx, UnitSystem.SI, UnitSystem.scaled)
dy = scaling.distance.convert(dy, UnitSystem.SI, UnitSystem.scaled)
print(f"{dx = :.2e} scaled")
print(f"{dy = :.2e} scaled")
# convert Fx and Fy to scaled units
Fx = scaling.water_vapor_flux.convert(Fx.values, UnitSystem.natural, UnitSystem.scaled)
Fy = scaling.water_vapor_flux.convert(Fy.values, UnitSystem.natural, UnitSystem.scaled)
# convert E to scaled units
# Do this for both the total E and the local regionally clipped E]
#print('pre-scaled',ds['Evap'])
E_total = scaling.evaporation.convert(ds["Evap"].values, UnitSystem.natural, UnitSystem.scaled)
E_local = scaling.evaporation.convert(ds["Evap_local"].values, UnitSystem.natural, UnitSystem.scaled)
L = 2.00e+06 m dx = 2.73e+04 m H = 2.00e+06 m dy = 2.74e+04 m dx = 1.37e-02 scaled dy = 1.37e-02 scaled
Plot the scaled E array
In [16]:
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS)
ax.add_feature(cartopy.feature.RIVERS)
ax.add_feature(cartopy.feature.OCEAN)
#Plot data
Ea = np.transpose(np.average(E_local,axis=2))
collection = plt.pcolormesh(ds['Evap'].lon, ds['Evap'].lat,Ea)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='black', alpha=0.5, linestyle='dotted')
gl.top_labels = False
gl.right_labels = False
#ax.set_title()
#plt.savefig()
plt.show()
plt.clf()
plt.close()
Run recycling model for each timestep
- Create recycling output array based on the shape of one of the surface input files: evap
- Translate evap and fluxes to secondary grid
- Calculate modeled precipitation
- Plot scaled input variables (evap and fluxes)
- Run through each timestep in the input files and calculate recycling ratio at each timestep across domain
- Plot rho and convergence metric for each timestep
In [17]:
import matplotlib.pyplot as plt
import logging
logging.basicConfig()
logging.getLogger("bulk_recycling_model").setLevel(logging.INFO)
from bulk_recycling_model import plotting
from bulk_recycling_model.main import run
#Make the rho array the same shape as the total E - will clip the external points at the end
rho_ar = np.empty((np.shape(E_total)[0]-1,np.shape(E_total)[1]-1,np.shape(E_total)[2]))
#Entering preprocessing and time step loop
#Run model and plot
for i,time in enumerate(ds.time):
# preprocess E onto the secondary grid
Ei_total = ED_preprocess.prepare_E(E_total[:,:,i])
Ei_local = ED_preprocess.prepare_E(E_local[:,:,i])
# preprocess water vapor fluxes onto the secondary grid
Fxi_left = preprocess.prepare_Fx_left(Fx[:,:,i])
Fxi_right = preprocess.prepare_Fx_right(Fx[:,:,i])
Fyi_bottom = preprocess.prepare_Fy_bottom(Fy[:,:,i])
Fyi_top = preprocess.prepare_Fy_top(Fy[:,:,i])
# compute P
Pi = preprocess.calculate_precipitation(Fxi_left, Fxi_right, Fyi_bottom, Fyi_top, Ei_total, dx, dy)
# # Create a quiver plot
# fig, ax = plt.subplots()
# U,V = plotting.build_uv_fluxes(Fxi_left, Fxi_right, Fyi_bottom, Fyi_top)
# X, Y = np.meshgrid(lon_axis.half_step, lon_axis.half_step, indexing="ij")
# ax.quiver(X[::2, ::2],Y[::2, ::2],U[::2, ::2],V[::2, ::2])
# fig.suptitle("Water Vapor Fluxes on cell edges")
# Create a quiver plot
fig, ax = plt.subplots()
collection = plotting.pcolormesh(ax, Ei_local, lon_axis, lat_axis, alpha=0.5)
fig.colorbar(collection, label="E (scaled)")
plotting.quiver(ax, Fxi_left, Fxi_right, Fyi_bottom, Fyi_top, lon_axis, lat_axis)
fig.suptitle("Evaporation + Water Vapor Fluxes on cell edges")
# Create a precip plot
fig, ax = plt.subplots()
collection = plotting.pcolormesh(ax, Pi, lon_axis, lat_axis, alpha=0.5)
fig.colorbar(collection, label="P (calculated)")
plotting.quiver(ax, Fxi_left, Fxi_right, Fyi_bottom, Fyi_top, lon_axis, lat_axis)
fig.suptitle(" ")
# Run the model
status = run(
Fxi_left,
Fxi_right,
Fyi_bottom,
Fyi_top,
Ei_local,
Pi,
dx,
dy,
R=0.2,
R_1=0.2,
max_iter=300,
tol=1e-3,
)
#Print timestep and status (converged or not) and add rho to recycling ration array
print(i,time.values)
print(status['k'])
rho_ar[:,:,i] = status["rho"]
# plot each timestep
fig, ax = plt.subplots()
cmap=plt.cm.viridis
cmap.set_extremes(under='red', over='orange')
collection = plotting.pcolormesh(ax, status["rho"], lon_axis, lat_axis,
vmin=0.0, vmax=1,
cmap=cmap)
collection_e = plotting.pcolormesh(ax, Ei_local, lon_axis, lat_axis, alpha=0.2)
fig.colorbar(collection,extend='both')
fig.suptitle(str(time.values)+" $\\rho$")
#plt.savefig(datap+"rho_"+str(time.values)+".png")
plt.show()
plt.close()
# plot the convergence
deltas = status["deltas"]
fig, ax = plt.subplots()
ax.plot(deltas)
ax.set_title("Convergence")
ax.set_xlabel("Iteration")
plt.show()
plt.close()
WARNING:bulk_recycling_model.main:Did not converge in 300 iterations
0 1994-01-01T00:00:00.000000000 300
WARNING:bulk_recycling_model.main:Did not converge in 300 iterations
1 1994-02-01T00:00:00.000000000 300
INFO:bulk_recycling_model.main:Converged in 211 iterations and 0:00:02.697145
2 1994-03-01T00:00:00.000000000 211
WARNING:bulk_recycling_model.main:Did not converge in 300 iterations
3 1994-04-01T00:00:00.000000000 300
INFO:bulk_recycling_model.main:Converged in 226 iterations and 0:00:02.803480
4 1994-05-01T00:00:00.000000000 226
WARNING:bulk_recycling_model.main:Did not converge in 300 iterations
5 1994-06-01T00:00:00.000000000 300
WARNING:bulk_recycling_model.main:Did not converge in 300 iterations
6 1994-07-01T00:00:00.000000000 300
INFO:bulk_recycling_model.main:Converged in 241 iterations and 0:00:03.180749
7 1994-08-01T00:00:00.000000000 241
WARNING:bulk_recycling_model.main:Did not converge in 300 iterations
8 1994-09-01T00:00:00.000000000 300
INFO:bulk_recycling_model.main:Converged in 196 iterations and 0:00:02.566546
9 1994-10-01T00:00:00.000000000 196
INFO:bulk_recycling_model.main:Converged in 260 iterations and 0:00:03.546760
10 1994-11-01T00:00:00.000000000 260
INFO:bulk_recycling_model.main:Converged in 256 iterations and 0:00:03.152958
11 1994-12-01T00:00:00.000000000 256
Create and save rho xarray file
- Create an xarray to store all of the calculated recycling ratios that is organised in an easy to plot/interpret format
- Count number of values in array over 1 - replace all of these with 1
- Count number of negative rho values - replace all of these with zero
- Save to file
In [18]:
lon_ar = np.linspace(start=ds.coords["lon"].min().values+lon_axis.step/2,
stop=ds.coords["lon"].max().values-lon_axis.step/2,
num=lon_axis.n_points-1)
lat_ar = np.linspace(start=ds.coords["lat"].min().values+lat_axis.step/2,
stop=ds.coords["lat"].max().values-lat_axis.step/2,
num=lat_axis.n_points-1)
rho_xarr = xr.Dataset(
data_vars=dict(rho=(["lon","lat","time"],rho_ar)),
coords=dict(
lon=(["lon"], lon_ar),
lat=(["lat"], lat_ar),
time=(["time"],ds.time.data)
),
attrs=dict(
description="Recycling ratio",
units="%",
),
)
rho_xarr = rho_xarr.transpose("time","lat","lon")
#rho_xarr = rho_xarr.sel(lat=slice(-9,3),lon=slice(18,30))
mask_lon = (rho_xarr.lon >= 18) & (rho_xarr.lon <= 30)
mask_lat = (rho_xarr.lat >= -9) & (rho_xarr.lat <= 3)
rho_xarr = rho_xarr.where(mask_lon & mask_lat, 0.0)
rho_xarr.to_netcdf(datao+"drc_rho_era5_"+str(YR)+".nc")
In [19]:
#Filtering out outliers for plotting
print('Number of rhos over 1: ', rho_xarr['rho'].where(rho_xarr['rho'].values>1.0).count().values)
print('Number of negative rhos: ', rho_xarr['rho'].where(rho_xarr['rho'].values<0.0).count().values)
rho_xarr = rho_xarr.where(rho_xarr['rho'].values<=1.0,1.0)
rho_xarr = rho_xarr.where(rho_xarr['rho'].values>0.0,0.0)
print('Number of rhos over 1: ', rho_xarr['rho'].where(rho_xarr['rho'].values>1.0).count().values)
print('Number of negative rhos: ', rho_xarr['rho'].where(rho_xarr['rho'].values<0.0).count().values)
end_all = timer.time()
length = end_all - start_all
print("Running the whole prep and recycling code took ", length, "seconds")
Number of rhos over 1: 0 Number of negative rhos: 32 Number of rhos over 1: 0 Number of negative rhos: 0 Running the whole prep and recycling code took 50.04628396034241 seconds
Plotting
Create seasonal arrays and plot these
In [20]:
mam_rho = rho_xarr['rho'].sel(time=rho_xarr.time.dt.month.isin([3,4,5]))
fig, ax = plt.subplots()
collection = mam_rho.mean("time").plot.contourf(vmin=0.0,vmax=0.6,levels=13,ax=ax,extend='max')
fig.suptitle("MAM "+str(YR)+" $\\rho$")
#plt.savefig(datao+"rho_MAM_"+str(YR)+".png")
plt.show()
son_rho = rho_xarr['rho'].sel(time=rho_xarr.time.dt.month.isin([9,10,11]))
fig, ax = plt.subplots()
collection = son_rho.mean("time").plot.contourf(vmin=0.0,vmax=0.6,levels=13,ax=ax,extend='max')
fig.suptitle("SON "+str(YR)+" $\\rho$")
#plt.savefig(datap+"rho_SON_"+str(YR)+".png")
plt.show()
jja_rho = rho_xarr['rho'].sel(time=rho_xarr.time.dt.month.isin([6,7,8]))
fig, ax = plt.subplots()
collection = jja_rho.mean("time").plot.contourf(vmin=0.0,vmax=0.6,levels=13,ax=ax,extend='max')
fig.suptitle("JJA "+str(YR)+" $\\rho$")
#plt.savefig(datap+"rho_JJA_"+str(YR)+".png")
plt.show()
djf_rho = rho_xarr['rho'].sel(time=rho_xarr.time.dt.month.isin([12,1,2]))
fig, ax = plt.subplots()
collection = djf_rho.mean("time").plot.contourf(vmin=0.0,vmax=0.6,levels=13,ax=ax,extend='max')
fig.suptitle("DJF "+str(YR)+" $\\rho$")
#plt.savefig(datap+"rho_DJF_"+str(YR)+".png")
plt.show()
plt.clf()
plt.close()
rho_xarr.close()